Dynamical Response of Networks under External Perturbations: Exact Results 
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We introduce and solve a general model of dynamic response under external perturbations. This 
model captures a wide range of systems out of equilibrium including Ising models of physical sys- 
tems, social opinions, and population genetics. The distribution of states under perturbation and 
relaxation process reflects two regimes — one driven by the external perturbation, and one driven by 
internal ordering. These regimes parallel the disordered and ordered regimes of equilibrium physical 
systems driven by thermal perturbations but here are shown to be relevant for non-thermal and 
non-equilibrium external influences on complex biological and social systems. We extend our results 
to a wide range of network topologies by introducing an effective strength of external perturbation 
by analytic mean-field approximation. Simulations show this generalization is remarkably accurate 
for many topologies of current interest in describing real systems. 
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Networks have become a standard model for a wealth 
of complex systems, from physics to social sciences to bi- 
ology [1, 2j. A large body of work has investigated topo- 
logical properties H, |j, @ . The raison d'etre, though, 
of complex network studies is to understand the rela- 
tionship between structure and dynamics - from disease 
spreading and social influence @, 0, [1, H, H3| to search 11 1. 
Yet, dynamic response of networks under external pertur- 
bations has been less thoroughly investigated [!, 121. In 
this paper we consider a simple dynamical process as a 
general framework for the dynamic response of a network 
to an external environment. The environment is initially 
treated as part of the network and then generalized as an 
external system. 

We obtain complete and exact results for the simplest 
case of fully connected networks and find a nontrivial dy- 
namic behavior that can be divided into two regimes. For 
large perturbations the environmental influence extends 
into the system with a distribution which, in the ther- 
modynamic limit, becomes a Gaussian around a value 
that reflects a balance between the external perturba- 
tions driving the system in different directions. For small 
perturbations the distribution of states has peaks at the 
two ordered states. Order arises from interactions within 
the system, and power law tails result from the external 
perturbation away from these ordered states. The bound- 
ary between these regimes is characterized by a uniform 
distribution where all states are equally likely. The time 
scale of equilibration is small for large perturbations and 
diverges inversely as the strength of the perturbation for 
small perturbations. This characterizes the switching 
time behavior of the two ordered states. We generalize 
the exact results to networks of different topologies using 
a mean field treatment. Simulations show that this gen- 
eralization, which involves renormalizing the constants 
in the distributions, is very accurate. Our results reveal 
and generalize key features of relaxation and dynamic 
response of models of a wide range of physical systems 



in the Ising universality class, electoral and contagion 
models of social systems, and the Wright-Fisher model 
of evolution in population biology. 

Specifically, we consider networks with N + No + N\ 
nodes. Each node has an internal state which can take 
only the values or 1. We let the Nq nodes be frozen 
in state 0, and N\ in state 1, and the remaining nodes 
change by adopting the state of a connected node. At 
each time step a random free node is selected; with prob- 
ability 1 — p the node copies the state of one of its con- 
nected neighbors, and with probability p the state re- 
mains unchanged. The frozen nodes can be interpreted 
as external perturbations to the subnetwork of free nodes. 
Analytically extending No and Ni to be smaller than 1 
enables modeling the case of weak coupling. This model 
generalizes our previous efforts to derive exact results of 
network dynamics [l3[ (see alsoflij). 

This system is similar to the Ising model, where 
Ni + No by explicitly representing the impact of thermal 
perturbations play the roles of the temperature T, and 
Ni—Nq acts as an external magnetic field h. Our dynam- 
ics are equivalent to Glauber dynamics [15| for weak fields 
and high temperatures, where the Ising model parame- 
ters are J/kT -> l/(z + N + Ni) and h/J -> (Ni-N ), 
where z is the number of nearest neighbors and J the 
nearest-neighbor interaction strength. For low tempera- 
tures our model is an alternative dynamics that also cap- 
tures the key kinetic properties of this system. Relevant 
network structures include crystalline 3-D lattices and 
random networks for amorphous spin-glasses; fully con- 
nected networks correspond to long range interactions or 
the mean field approximation. Despite the relevance to 
the extensively studied Ising model, we are not aware of 
any other exact solution of the response dynamics of a 
fully connected system or explicit representation of ther- 
mal or other perturbation for dynamic response. Specific 
results are available only for zero temperature dynamics 
in one-dimensional or mean field systems. (l6l . [l7j 
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Our system can also model an election with two can- 
didates jlH [l9| where some of the voters have a fixed 
opinion while the rest change their intention according 
to the opinion of others. Another application is to epi- 
demics that spread upon contact between infected nodes 
(e.g., individuals or computers). Finally, the model can 
represent an evolving population of sexually reproducing 
(haploid) organisms where the internal state represents 
one of two alleles of a gene (20]. Taking p — 1/2, the 
update of a node mimics the mating of two individuals, 
with one parent being replaced by the offspring, which 
can receive the allele of either the mother or the father 
with 50% probability. Since a free node can also copy the 
state of a frozen node, the ratios Nq/(N + No + Ni — 1) 
and Nx/(N + No + Ni — 1) give the mutation rates. 

For a fully connected network the nodes are indistin- 
guishable and the state of the network is fully specified 
by the number of nodes with internal state 1 [13| . There- 
fore, there are only N + 1 global states, which we denote 
k = 0, 1, N. The state has k free nodes in state 
1 and N — k free nodes in state 0. If Pt(m) is the proba- 
bility of finding the network in the state o~ m at the time 
t, then P t+ i(m) can depend only on P t (m), Pt(m + 1) 
and Pt(m — 1). The probabilities Pt(m) define a vector 
of N + 1 components P*. In terms of Pt the dynamics is 
described by the equation 



P t+ i = UP t = 1 



(I-P) 



A P, 



N(N + N + Ni - 1) 

where the evolution matrix U, and also the auxiliary ma- 
trix A, is tri-diagonal. The non-zero elements of A are 
independent of p and are given by 

A m . m = 2m(N -m)+Ni(N- m) + N m 
A m , m +i = -(m + 1)(N + N - m - 1) 
A m ,m-x = -(N - m + l)(Ni + m - 1). 

The transition probability from state o~m to o~l after a 
time t can be written as 



P(L,t;M,0) 
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n r 
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(1) 



where a r L and b r u are the components of the right and 
left r-th eigenvectors of the evolution matrix, a r and b r , 
with r r = b r a r . Thus, the dynamical problem has been 
reduced to finding the right and left eigenvectors and the 
eigenvalues of A. 

It is easy to check by inspection of small matrices that 
the eigenvalues /i r of A are given by 



r(r- I + Nq + N^. 



This implies < p < X r < 1 , where X r are the eigenvalues 
of U. Because of Eq.([T]), the unit eigenvalues completely 
determine the asymptotic behavior of the system. 

The eigensystem Aa r = p r & r leads to the following 
recursion relation for the coefficients a rm 



m+l 



A, 



j=m-l 



j-iif (Xy 
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with a r! jv+i = Gtr,-i = 0. To solve this equation we mul- 
tiply the whole expression by x m , sum over m and de- 
fine the generating function p r (x) — 53m=o a rmX m ■ 

The 

recursion relation then yields the following differential 
equation for p r 

x{l - x)p'l + [(1 - N - N ) - (1 + Ni - N)x]p' r + 
[JVJVi - Mr/(1 - x)]Pr = 0. 

(3) 

To understand the asymptotic behavior of the system 
(/i r = 0) we have to consider two cases: 

(a) If N = Ni = then fi r = leads to r = or r = 1 
[jjj ]. In this case the differential equation simplifies to 
xp'l + (1 — N)p' r = 0, whose two independent solutions 
are po(x) = 1 and pi(x) — x , corresponding to the 
all-nodes-0 or all-nodes-1 states respectively. 

(b) H N ,Ni ^ then fx r = implies r = 0. In this 
case equation (J3|) is that of a hypergeometric function F 
and we find po(aO = F(—N, Ni,l — N—No,x), which is a 
finite polynomial with known coefficients ao m ■ Normaliz- 
ing this eigenvector, we obtain the probability of finding 
the network in state o~ m at large times: 

[N — my. m\ 

where A = A(N, No, Ni) is a normalization. Because 
of the frozen nodes, the dynamics will never stabilize in 
any state, but will always move from one state to another, 
with mean occupation number m = NNi/ (No+Ni). The 
surprising feature of this solution is that for TVq = = 1 
we obtain p{m) = 1/(N + 1), for all values of N. Thus 
all macroscopic states are equally likely and the system 
executes a random walk through the state space. 

The dynamics at long times is dominated by the sec- 
ond largest eigenvector, with eigenvalue Ai. For large 
networks ~ e -t / T where 



N(N + N + Nx-1) 

(l-p^No + m) 



(•5) 



We obtain a complete description of the dynamics by 
deriving all eigenvectors with /x r ^ 0. The differential 
equation for p r {x) yields 

/ \ _ Fjl-r-NaA-r-N-No-N^l-N-No^) / r n 
Pr\x) — h _ a ,)r-l+JV +JV 1 • {V) 

Expanding the numerator and denominator in Taylor 
series gives the coefficients a rm . Although they can 
easily be written down explicitly, we do not do so 
here. Similarly, defining the generating function q r (x) = 

Sm=i^./Vi brm.x m w & obtain a differential equation for q r 
whose solution is 



q r (x) 



"2 F(l-r-N 1 .l-r-N-N n -N 1 ,l-N-N 1 .x) 
(l-x)H- 1 ' 



If A^q = Ni — this solution is not valid for r = or 
r = 1, since the matrix A T becomes singular. In this 
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FIG. 1: (color online) Asymptotic probability distribution for 
a network with N = 100 and several values of iVo and N\. 



case the two left eigenvectors are given by &o,m = 1 and 
bi.m — N — 2m. For other cases the solution is obtained 
from the power series expansion of q r (x). Equations © 
and complete the solution of the problem. 

In the thermodynamic limit N — > oo we can define 
continuous variables x — m/N, no = Nq/N and n\ — 
Ni/N and approximate the asymptotic distribution by 
a Gaussian p(x) — po exp [— (x — xq) 2 /2S 2 ] with xq = 



Til I (no + m), po = 1/ v 2ttS 2 and 



S = 



n ni(l + n + n\) 



N(n + niY 



1/2 



(8) 



In the limit where no, ni » 1 the width depends only 
on the ratio a — no/n\ and is given by y/a/N /(l + a). 

The problem we just solved can be generalized to treat 
an external reservoir weakly coupled to the network of N 
nodes. We note that the differential equations for the 
generating functions p r (x) and q r (x) remain well defined 
for real No and N\. The solutions for the generating 
functions remain the same, except that factorials must be 
replaced by gamma functions. Since Nq/(N+Nq+Nx—1) 
and N\/(N + No + Ni — 1) represent the probabilities 
that a free node copies one of the frozen nodes, small 
values of No and N\ can be interpreted as representing a 
weak connection between the free nodes and an external 
system containing the frozen nodes. The external system 
can be thought of as a reservoir that affects the network 
but is not affected by it. Alternatively, we can suppose 
that there is a single node fixed at that is on for only 
a fraction No of the time and off for the fraction 1 — No, 
and similarly for a single node fixed at 1. 

Figure 1 shows examples of the distribution p(m) for 
a network with N = 100 and various values of Nq and 
N\. Numerical simulations displaying similar results are 
described in [2lj . 

Figure 2 shows an example of the time evolution of the 
probability density for a fully connected network com- 
pared to numerical simulations. The evolution from the 



FIG. 2: (color online) Time evolution of the probability distri- 
bution P t for a network with N = 100 and N = Nx = 5. The 
initial state is Poi = 0.5(<5i,2o + <5i,so). The histograms show 
the average over 50,000 actual realizations of the dynamics 
and the solid (red) line shows the analytical result. 



initial to the asymptotic time- independent distribution is 
the analog of an equilibration process promoted by the 
external system. 

For small values of Nq and Ni (<< 1/lnAT), we can 
obtain a simplified expression for p(m): 



p{m) 



N x Np 
N + N ± 



1-JVxlnJV 1-JVqIhJV 



,1-N! 



(N 



(9) 



Thus p(m) displays a power law behavior on both ends of 
the curve: 1/m for m close to and 1 / (N—m) for m close 
to N (sec, for instance, the curve with No = N\ = 0.5 
in Fig. 1). Since the relaxation time r is proportional 
to 1/(Nq + Ni), the equilibration process becomes very 
slow in this limit. 

For networks with different topologies the effect of the 
frozen nodes is amplified. To see this we note that the 
probability that a free node copies a frozen node is P\ — 
(No + N 1 )/(No + Ni + la) where h is the degree of the 
node. For fully connected networks ki = N — 1 and 
Pi = Pfc- For general networks an average value P av can 
be calculated by replacing fcj by the average degree k av . 
We can then define effective numbers of frozen nodes, 
Noef and Ni e f , as being the values of No and N% in Ppc 
for which P av = Pfc- This leads to 



Noef = fN , 



Nut - fNi 



(10) 



where / = (N — 1) / k av . Corrections involving higher mo- 
ments can be obtained by integrating Pj with the degree 
distribution and expanding around k av . 

Figure 3 shows examples of the equilibrium distri- 
bution for four different networks with N = 100 and 
No = Ni = 5. Panel (a) shows a random network 
with connection probability 0.3 (N av — 30, / = 3.3). 
The theoretical result was obtained with Eq. 0$ with 
Noef = N\ e f = 17. For a scale-free network (panel (b)) 
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FIG. 3: (color online) Asymptotic probability distribution for 
networks with different topologies. In all cases N = 100, No = 
Ni = 5, t = 10, 000, and the number of realizations is 50, 000. 
The theoretical (red) curve is drawn with effective numbers 
of frozen nodes Noef = /JVn and -/Vie/ = fNi: (a) random 
network No e f = N\ e f = 17; (b) scale-free No e f ~ Ni e f = 82; 
(c) regular 2-D lattice iVn e / = Ni e f = 140; (d) small world 
network N 0e f = N lef = 140. 



grown from an initial cluster of 6 nodes adding nodes with 
3 connections each following the preferential attachment 
rule [l|, / = 99/6 and the effective values of Nq and Nx 
are approximately 82. Panel (c) shows the probability 
distribution for a 2-D regular lattice with 10 x 10 nodes 
for which / = 99/3.6 ss 28. Finally, panel (d) shows a 
small world version of the regular lattice [l], where 30 
connections were randomly re-connected, creating short- 
cuts between otherwise distant nodes. These results show 
that the mean field generalization is accurate for many 
network topologies. Still, extreme cases such as a star 
network should be different and this is confirmed by sim- 
ulations and preliminary analytic results. The relaxation 



time ([5]), in units of network size and for p — 0, becomes 
t/N = (k av + N n + N 1 )/(N + iVi). It increases linearly 
with N for fully connected or random networks, but is 
independent of N for regular and scale-free topologies. 

Our results have important implications for real sys- 
tems. In the social sciences they show the importance 
of opinion makers in stabilizing the outcome of elections: 
weak external influences result in an arbitrary but seem- 
ingly strong opinion that can switch at random (see also 
|22|). due to the arbitrary choice of the ordered state in 
the weak perturbation regime. Thus, an elected candi- 
date winning a landslide election may have no solid sup- 
port. The slow dynamics can play a crucial role, since 
the time to switching might occur only after the election 
day, especially if the number of voters is large. Stronger 
external influences, counter intuitively, reduce the rela- 
tion time giving rise to improved internal equilibration. 
In theoretical biology our results are equivalent to the ex- 
act dynamical solution of the Wright-Fisher model [2(| 
for arbitrary population sizes and mutation rates. Our 
equations give not only the asymptotic equilibrium dis- 
tribution of alleles (see }20| for approximate expressions) , 
but also its time evolution in two regimes, one where mu- 
tations have difficulty overcoming an existing dominant 
allele (the low perturbation regime) and one where ran- 
dom mutations dominate (the high perturbation regime) . 
Again this is crucial information, since the equilibration 
time can be extremely long for the typically small muta- 
tion rates observed in nature. 

Finally, we emphasize that exact dynamical solutions, 
describing systems out of equilibrium, are rare and 
important for the study of many statistical properties 
that can be described in future publications. 
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